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Abstract 

We consider lattice self-avoiding walks and discuss the dynamic critical behavior 
of two dynamics that use local and bilocal moves and generalize the usual reptation 
dynamics. We determine the integrated and exponential autocorrelation times for 
several observables, perform a dynamic finite-size scaling study of the autocorrelation 
functions, and compute the associated dynamic critical exponents z. For the variables 
that describe the size of the walks, in the absence of interactions we find z ~ 2.2 in 
two dimensions and z ~ 2.1 in three dimensions. At the ^-point in two dimensions we 
have z ~ 2.3. 



1 Introduction 



The lattice self-avoiding walk (SAW) is a well-known model for the critical behavior of a 
homopolymer in a solvent and it has also been extensively used in the study of several 
properties of heteropolymers The earliest simulations either used a local dynamics 

in which, at each step, a small part of the walk (usually 2-4 consecutive beads) was modified, 
or the so-called reptation dynamics 1^-0. All these algorithms are however nonergodic ||^,^ 



r0|-[12[| and only a limited fraction of the phase space is visited. Note that, contrary to some 



claims in the literature, the deviations are sizeable even for very short walks if one is interested 
in low-temperature properties, i.e. polymers in the compact phase or heteropolymers near 
the folding temperature ||I3 H15|] . For instance — see footnote 4 in Ref. |T6| — if one uses the 
Verdier-Stockmayer algorithm in two dimensions, one does not sample 3.2%, 1.4%, 5% 
of the most compact configurations for = 11, 13, 15. These ergodicity problems can be 



solved by using a different ensemble [p]^-pO||, chain-growth algorithms pl|-p4|, or nonlocal 
algorithms P5|-PD[]. However, in the presence of interactions, nonlocal algorithms become 
inefficient since nonlocal moves generate new walks with large energy differences and thus 
they will be rejected making the dynamics very slow. Chain-growth algorithms may work 
much better but, in order to make them efficient, one must sample a different probability 
distribution, and thus one may be worried by the introduced bias. 

In this paper we wish to discuss a family of dynamics that use bilocal moves and generalize 
the reptation dynamics: a bilocal move alters at the same time two disjoint small groups of 
consecutive sites of the walk that may be very far away. Since a small number of beads is 
changed at each step, these algorithms should be efficient in the presence of interactions, and 
thus they can be used in the study of the collapsed phase and of the folding of heteropolymers. 
Similar moves were introduced in Refs. [01| 



The ergodicity properties of these algorithms have been discussed in Ref. |[16|. Here, 
we will study the dynamic properties of two different implementations. The first one, the 
extended end-end reptation (EER) algorithm, is obtained by performing at the same time 
reptation moves and bilocal kink- kink moves. Such an algorithm is quite efficient. In the 
absence of interactions, the autocorrelation time for quantities that measure the size of 
the walk — for instance, the end-to-end distance or the radius of gyration — scales as A^^ 
with z ~ 2.2 in two dimensions and z ^ 2.1 in three dimensions. The behavior of the 
energy (number of nearest-neighbor contacts among nonconsecutive links) is even faster, 
with 2; ^ 1.7 in both dimensions. We have also tested the behavior of the algorithm at the 
Q point in two dimensions. We find that the critical behavior is only marginally worse, with 
z ^ 2.3, both for metric quantities and for the energy. We have also studied a different 
version, the extended kink-end reptation (KER) algorithm, in which the reptation moves 
are replaced by kink-end moves, i.e. by moves in which a kink is cleaved and two additional 
links are attached at the end of the walk and viceversa. This version is much slower, with 
z ^ 2.9 for quantities that measure the polymer size. Clearly, reptation moves are essential 
to obtain a fast dynamics. 

The paper is organized as follows. In Sec. ^ we define the model and the observables 
whose critical behavior will be studied. In Sec. ^ we define the basic moves and the two 
dynamics. Specific implementation details are reported in the Appendix. In Sec. § we 
define the autocorrelation times, the dynamic critical exponents, and the methods we use to 
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compute them. In Sees. ^ and |^ we discuss the critical behavior of the EER and of the KER 
algorithms in the absence of interactions, while in Sec. ^ we discuss the behavior of the EER 
algorithm in two dimensions at the 9 point. Conclusions are presented in Sec. ||. 



2 Definitions 

In this paper we consider SAWs with fixed number of steps N and free endpoints on a 
hypercubic lattice U^. An A^-step SAW a; is a set of A^ + 1 lattice sites cuq, . . ., cutv, such that 
UJi and cjj+i are lattice nearest neighbors. By translation invariance we may fix ujq to be the 
origin. 

We consider several observables that measure of the size of an A^-step SAW: 

• The mean square end-to-end distance 

Rl = {uoN-uJof ■ (1) 

• The mean square radius of gyration 

1 ^ / 1 ^ \' 1 ^ 

Rl = y^U^i y^cjfc =— —S^iuOi-uo.f. (2) 

^ N + J 2(N +iy ^' ^ ' 

• The mean square monomer distance from an endpoint 



1 ^ 

m = T^— tXI^^^ (3) 

i=0 



R 



Moreover, for each SAW, we define the number of nearest-neighbor contacts S. It is defined 
as follows. For every i ^ j we define 





if \uji — uij\ 




otherwise. 



N-2 N 



- 1 n ;„„ W 

Then, 

^ - - E E (5) 

i=0 j=i+2 

Note that we do not include here the trivial contacts between consecutive walk sites. 
We give each walk u a weight W{(3) depending on the inverse temperature (3 = 1/kT, 

W{f3) ^ ^e-'^^ (6) 

where is the partition sum 

Zn = J2^-'', (7) 
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square lattice 







year 


method 




Ref. 




1987 


EE 


0.75 


Ref. 


43 


1988 


MC 


0.65 ±0.03 


Ref. 


44 


1990 


EE 


0.67 ±0.04 


Ref. 


45 


1990 


MC 


0.658 ±0.004 


Ref. 


46 1 1992 


EE 


0.657 ±0.016 


Ref. 


47 


1993 


MC 


0.658 ±0.004 


Ref. 


48 


1994 


EE 


0.660 ±0.005 


Ref. 


49 


1995 


MC 


0.665 ±0.002 


Ref. 


5q 


1996 


MC 


0.664, 0.666 


cubic lattice 






year 


method 




Ref. 


51] 1989 


MC 


0.274 ±0.006 


Ref. 


52 


1995 


MC 


0.2687 ±0.0005 


Ref. 


53 


1996 


MC 


0.2779 ±0.0041 


Ref. 


53 


1996 


MC 


0.2782 ±0.0070 


Ref. 


54 


1996 


MC 


0.276 ±0.006 



Table 1: Estimates of I3e on the square and cubic lattices. EE stands for exact enumeration, 
MC for Monte Carlo. 



and the sum is extended over all A^-step SAWs. The mean values {RI)n, {RI)n, and {R^) 



N 



have the asymptotic behavior 

{RDn, {RI)n, {RDn ~ N"" (8) 

as N oo, where i/ is a critical exponent, which is independent of the microscopic details 
of the lattice model but depends on /?. For P < Pe (good-solvent regime) the SAW is swollen 



and (see Ref. for a review of estimates of u in three dimensions) 

d = 2. 




(9) 

d = 3 (Ref. [36]), 



while for P > the SAW is compact with 
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(A) 



(B) 



(C) 



Figure 1: All one-bead moves: (A) One-bead flip. (B) 90° end-bond rotation. (C) 180'^ 
end-bond rotation. 



For the very specific value P = Pe (6'-point) 

4 

X log 



d = 2 (Ref. [37]); 



d = 3 (Refs. []I|,§|-|T 



The critical value Pe depends on the model, on the lattice type, and on all microscopic 
details. For the model we consider here on the square and on the cubic lattice, the best 
estimates are reported in Table IT]. 



3 Algorithms 

For the simulation of weakly interacting walks, i.e. for /5 < there exist powerful nonlocal 



algorithms [^^|30|. However, these algorithms cannot be used in confined geometries — they 
are not ergodic — and are very inefficient in the presence of strong interactions. Indeed, in 
these conditions nonlocal moves are rarely accepted. In this paper we consider two algorithms 
that use local and bilocal moves . A local move is one that alters only a few consecutive 
beads of the SAW, leaving the other sites unchanged. A bilocal move is instead one that 
alters two disjoint small groups of consecutive sites of the walk; these two groups may in 
general be very far from each other. 

In our study we consider two types of local moves (see Fig. [|): 

[LO] One-bead flips in which one internal bead (i.e. uj{i) , 1 < i < N — 1) only is moved. 
[LI] End-bond rotations in which the last step of the walk is rotated. 
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C D 



C D 



A 



B 



A 



B 



Figure 2: The kink-transport move. A kink has been cleaved from AB and attached at CD. 
Note that the new kink is permitted to occupy one or both of the sites abandoned by the 
old kink. 




Figure 3: The kink-end reptation ( — >) and end-kink reptation {< — ) moves. In ( — >), a 
kink has been cleaved from AB and two new steps have been attached at the end marked X. 
Note that the new end steps are permitted to occupy one or both of the sites abandoned by 
the kink. 



We also introduce several types of bilocal moves: 

[B22] Kink-transport moves in which a kink is cleaved from the walk and attached at a pair 
of neighboring sites somewhere else along the walk (see Fig. note that the new kink 
is allowed to occupy one or both of the sites abandoned by the old kink. 

[BKE] Kink-end and end- kink reptation moves (see Fig. |^). In the kink-end reptation move 
a kink is deleted at one location along the walk and two new bonds are appended in 
arbitrary directions at the free endpoint of the walk. Viceversa, an end-kink reptation 
move consists in deleting two bonds from the end of the walk and in inserting a kink, 
in arbitrary orientation, at some location along the walk. 

[BEE] Reptation move (see Fig. ^ in which one bond is deleted from one end of the walk 
and a new bond is appended in arbitrary direction at the other end. 
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r 

I 



Figure 4: The reptation move. The head of the walk is indicated by X. The dashed hnes 
indicate the proposed new step and the abandoned old step. 



Using these elementary moves we introduce three different updates that leave invariant 
the Gibbs measure (^. They are described in detail in the Appendix. We will then consider 
two algorithms: 

(a) Extended end-end reptation (EER) algorithm; 

(b) Extended kink-end reptation (KER) algorithm. 

The EER algorithm consists in combining with non-zero probability the reptation move 
and the kink-kink local/bilocal move, see Appendices [A.l| and |A.2| for the implementation 



details. More precisely, with probability p one performs a reptation move, and with probabil- 
ity 1 — p a kink-kink local/bilocal move. The KER algorithm works analogously: instead of 
the reptation moves we use the kink-end/end-kink BKE moves. Both algorithms are known 
to be ergodic in two dimensions, while for d = ?> ergodicity has been proved only for the KER 



algorithm [T^. For the EER algorithm ergodicity is still an open problem. The probability 
p is not fixed, and can be tuned to obtain the fastest dynamics. However, it should not have 
any influence on the dynamic critical behavior, i.e. on the dynamic critical exponents. For 
this reason we have not performed a systematic study of the dependence of the numerical 
results on p and we have simply set p = 0.50 in most of our simulations. 

In order to have a fast and efficient implementation, it is important that local and bilocal 
moves are performed in a CPU time of order one, i.e. constant as N oo. This can be 
obtained only with a careful choice of the data structures used to store the walk coordinates. 
For local and end-end reptation moves it is sufficient to store the walk coordinates as a 
circular list in which the coordinates of the beads are stored sequentially. However, such a 
data structure is not convenient for B22 and BKE moves, since in this case insertion and 



deletion of a single point requires a time of order A^. The problem is solved |^5[ by storing 



the coordinates in a contiguously allocated doubly-linked linear list. The coordinates are 
stored contiguously but not in any particular order. In order to find the bead that follows 
and precedes a given one, one keeps two arrays of pointers that give the location in the 
coordinate list of the preceding and of the successive walk bead. With this type of data 
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structure, insertions and deletions take a time of order one. An efficient implementation 
requires also the ability to perform the self-avoidance check in a time of order one. This 
may be obtained by using a bit table — however, this is only possible for short walks — or a 
hash table. In this second case, one must be careful to use a hashing method that allows to 
insert and delete a single point in a time of order one. A very efficient method is described 



in Ref. |56| 



4 Dynamical behavior 



In order to determine the efficiency of the dynamics, we study the critical behavior of the 
autocorrelation times. Quite generally, given an observable A, we define the unnormalized 
autocorrelation function 

CAA{t) = {AsAs+t) - {A)\ (12) 



and its normalized counterpart 



PAA = CAA{t)/CAA{0), 



(13) 



where t is the dynamics "time." Typically, paa(^) decays exponentially, i.e. pAA{t) ~ e l*'/"^, 
for large \t\. Thus, we can define the ex])onentm/ autocorrelation time 



Texp,A = lim sup 



\t\ 



-\n\pAAit)\ 



(14) 



and 



Texp — SUpTexp,A- 

A 



(15) 



Thus, Texp is the relaxation time of the slowest mode in the system. 
We also define the integrated autocorrelation time 



Tint, A 



^ +00 ^ +CXD 



Paa{s) 



(16) 



s=l 



which controls the statistical error in Monte Carlo measurements of {A). The factor 1/2 is 
inserted so that rint,A ~ Te^p,A if PAAif) = e'l^l/^^-p-* with Texp.A > 1- 

In order to estimate t^^p^a, we will proceed as follows. If pAA{t) is the normalized auto- 
correlation function estimated from the data, we define an effective exponent 



exp, A 



{t;s) 



In 



pAAjt) 
pAAit+ S] 



;i7) 



where s is some fixed number. This quantity should become independent of t as t — ^ oo. 
In practice, we look for a region in which the estimates are sufficiently stable with t and 
then take the value of fexp,A(^; s) in this region as our estimate of Tgxp.A- In order to apply 
the method we should also choose the parameter s. It should be neither too small, i.e. 
s <^ Texp,A, otherwise Paa(^) ~ Paa(^ + s), nor too large, i.e. s ^ Tcxp,A5 otherwise the 
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error on pAA{t + s) is large. In our study we have always fixed s self-consistently, by taking 
s ~ Texp,A/k with k ^ 10-20. 

In order to estimate Tmt.A, "we use the self-consistent windowing method proposed in 
Ref. [0]. We define our estimate as 



M 

i ■ 

Tint, A 



1 ^ 



2 

t=-M 

where M is the smallest integer that satisfies M > cfint,A and c is a fixed constant. The 
variance of fi^t,A is given by 

, 2(2M+1) 2 . X 

Var Tint.A ~ r-^nt,A^ 19 

n 

where n is the number of Monte Carlo iterations and we have made the approximation 
Tint, A <^ M <^ The method provides quite accurate estimates of fi^t,A as long as c is 
chosen so that M is a few times r^xp^A- There are cases in which rint,A 'rexp,A, so that 
the previous condition is difficult to satisfy. In this case, a successful "ad hoc" recipe was 



proposed by Li et al. fS^] for the pivot algorithm. They noted that, when A is one of the 
radii, pAAif) ~ ^It^ in the intermediate region rint.A ^ ^ ^ Texp,A, with q ^ 1-1.2. Thus, 
they extrapolated Paa(^) proportionally to l/t for t > M, (i.e. set paa(^) = MpAA{M)/t for 
t > M), and then used this expression to compute the contribution to rint,A from the region 
M < t < Texp^A- This gives the modified estimate 

. +M 

fint,A = - 5^ PAAit) + MpAA{M) In (^) . (20) 

t=-M 

In general, the autocorrelation times diverge as oo. We can thus define two dynamic 

critical exponents -Zexp.A and ^int,A by 

T-int,A ~ AT^'"*'^, 

rexp,A ~ iV^-p.-^. (21) 

It is important to stress that in general Zint,A 7^ 2;exp,A and that different observables may 
have different critical exponents. 

The dynamic critical exponents can also be obtained by requiring the autocorrelation 
function to obey a dynamic scaling law of the form [p^,|58|,p9| 



PAA{t\N) ^ \t\-''FAA {tN-') ^ N-'^'GAAitN-") , (22) 

valid in the limit — 00, \t\ 00, with tN~^ fixed. Here, a, 6 > are dynamic critical 
exponents and Faa and Gaa are suitable scaling functions. If Faa{x) is continuous, strictly 
positive, and rapidly decaying — e.g. exponentially — as |x| 00, then, it is not hard to see 
that 

2:cxp,A = (23) 
-2int,A = (l-a)6 (if a < 1) . (24) 
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80000 
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120000 



160000 



Figure 5: Effective exponent fexp,R2(t;s) for the EER algorithm. Here d = 2, N 



s = 2000. The horizontal line corresponds to our final estimate Texp,ij2 
dashed lines to ± one error bar (900). 



300, 



43600 and the 



The exponents a and b can be determined by requiring the collapse onto a single curve 
of the autocorrelation functions corresponding to different values of A^. Then, using the 
previous formulae, one can determine Zexp,A and ^^mt.A- The advantage of this method is in 
its bypassing the problem of determining rexp,A and rint,A, but it is quite difficult to assess 
the errors, since the optimal values are determined visually. 



5 The EER dynamics in two and three dimensions 

We performed an extensive Monte Carlo simulation of noninteracting {f3 = 0) SAWs on 
the square lattice, d = 2, and on the cubic lattice, d = 3, using the EER algorithm. We 
set p = 0.5 and used the first version of the reptation move, see App. [A.2| . We considered 
only three values of N, N = 100, 300, 1000, but for each of them we collected a very large 
statistics, see Table |. We measured (i?g)jv, {RI)n, {RI,)n, and the energy £]\f, i.e. the 
number of nearest-neighbor contacts. We compared the static results of our simulations 



with those of Li et al. |57|, finding good agreement 
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d 


N 


iter 


1 exp,/c^ 


' exp,i?^ 


' exp,i?^;. 


1 expjfc 




100 


5.77 


10^*^ 


3230 ±20 


3062 ±24 


3190 ±16 


3078 ±40 


2 


300 


5.38 


10" 


43600 ±900 


47300 ± 1500 


43140 ±760 


35600 ±1260 




1000 


6.90 


10" 


710400 ±42000 


859600 ±46000 


745000 ±36000 


448000 ±30000 




100 


4.36 


10^" 


2980 ±40 


3020 ±30 


2968 ±30 


2840 ±60 


3 


300 


4.42 


10" 


32640 ± 840 


34700 ±900 


32980 ±840 


27480 ±800 




1000 


5.20 


10" 


370400 ±22000 


374800 ±24000 


368000 ±22000 


300800 ±16000 



Table 2: Exponential autocorrelation times for the EER algorithm in two and three dimen- 
sions, "iter" is the number of iterations. 



d 




'^exp,i?2 


'^exp,i?2 


'^cxp,i?^ 


'^exp,£ 






2.37±0.02 


2.46 ±0.02 


2.36±0.02 


2.18±0.02 


2 


B 


0.0610 ±0.0048 


0.0376 ±0.0032 


0.0586 ±0.0030 


0.133±0.013 






1.18 


1.03 


0.52 


2.33 






2.10±0.03 


2.11 ±0.04 


2.10±0.03 


2.04 ±0.03 


3 


B 


0.196 ±0.032 


0.180 ±0.028 


0.190±0.030 


0.242 ±0.028 






1.13 


1.90 


1.31 


1.57 



Table 3: Dynamic exponent Zexp for the EER algorithm in two and three dimensions, obtained 
by fitting r^xp = -BA^^^p. The number of degrees of freedom of the fit is 1. 



First of all, we measured the exponential autocorrelation times. As an example, in Fig. 
D we report the effective exponent Texp,_R2(t; s) for = 300, d = 2, s = 2000. It increases 
with t and becomes constant within error bars for t > 106000. We have taken our estimate 
at t ~ 130000 ~ 3rexp,R2. The same analysis has been done for all observables and all values 
of A^. The results are reported in Table 0. 

We have determined the dynamic exponent z^xp.A by fitting the results of Table ^ with 
the Ansatz 

rexp,A = BaN'^--'^. (25) 

The results are reported in Table |^. Since we only have three values of A^, we cannot study 
the effect of corrections to scaling and thus we cannot determine the systematic error on our 
results. However, an indication can be obtained by comparing the results for the three radii. 
Indeed, one expects these three observables to have the same dynamic exponent. In two 
dimensions, 2;exp,ij2 differs significantly from the estimates for the other radii, indicating the 
presence of large corrections to scaling. Such conclusion is confirmed by the scaling analysis 
that will be reported below: as we shall see, the results of Table |^ in two dimensions are 
effective exponents that are expected to decrease as larger values of A^ are included in the fit. 
Of course, it is possible that also Zexp,^' is affected by large corrections to scaling. However, 
the scaling analysis confirms the result for Zexp,^'; thus showing that the scaling corrections 
are important for the radii only. 

In three dimensions, there is no evidence of large scaling corrections. Indeed, the expo- 
nents Zexp for the three radii agree within error bars. The dynamic exponent for the energy 
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is lower than that for the radii. However, the difference is small and it is not clear if it is 
real or just the result of neglected corrections to scaling. 

Then, we have computed the integrated autocorrelation times by using the self-consistent 
windowing method with c = 15, cf. Eq. ([T8|). The results are reported in Table |[ 

For the radii, Texp,R2 ~ 1-2 rint,i?2 and thus the choice c = 15 should be rather conservative. 
For the energy, c = 15 corresponds to 1.3rexp,£ ^ M < 4rexp,f in two dimensions and to 
2 Tcxp,£' ^ M < 5 Texp,e in three dimensions for our values of (M is the cutoff defined in Eq. 
([T8|)). Such values of M can give rise, at least for = 1000, to a systematic underestimate of 
rint,£. For our choice of c and for = 1000, peeiM) ^ (7.2 ± 1.5) • 10"^ and (4.9 ± 1.0) • 10"^ 
in two and three dimensions respectively. Therefore, assuming a pure exponential behavior 
for if: > M, the neglected contribution should be of order 320 ± 60 in c? = 2 and 150 ± 30 
in (i = 3. The correction is quite small, but larger than the quoted error bars: the correct 
errors are at least larger by a factor of 5 (resp. 2) in two (resp. three) dimensions. 

The results for rint,A have been fitted with the Ansatz 

rint,A = BaN'''^^^\ (26) 

The results are reported in Table ^. 

The error bars are purely statistical and, as in the case of Zexp, large systematic errors 
may be present. We expect the quantities that measure the walk size to have the same 
dynamic exponent Zinf Thus, both in two and three dimensions, the reported errors are 
largely underestimated. More conservative estimates are 

Zint,/?2 = 2.20 ± 0.03 d = 2, 

Zint,R^ = 2.07 ± 0.02 d = 3. (27) 

The fit for the energy has a very large x^- There are two reasons for this. First, the statistical 
errors are underestimated, as we already discussed. Second, there may be large — compared 
to the tiny statistical errors — corrections to scaling. For these reasons, the errors quoted 
in Table ^ for Zint,£ should not be taken seriously. A more realistic estimate is obtained by 
multiplying the errors by which gives an error of ±0.04 on the exponent in both two 

and three dimensions. Such error is more realistic and indeed is close to the error we quoted 
for the radii, cf. Eq. (p^). 



d 


N 




int,i?2 




Tmt,£ 




100 


3113.6±5.6 


2209.2 ±3.4 


2124.4±3.2 


861.64±0.82 


2 


300 


36000 ±72 


24250 ±40 


23764 ±38 


5596.4 ±4.4 




1000 


522880 ±3520 


337220 ±1820 


334600 ±1800 


39224 ±72 




100 


2729. 0±5. 2 


1761. 2±2. 8 


1780.0 ±2.8 


925.4±1.0 


3 


300 


27092 ±52 


16772 ±26 


17162 ±26 


6035.8 ±5.4 




1000 


332900 ±2100 


199160 ±960 


205100 ±1000 


42954 ±96 



Table 4: Integrated autocorrelation times for the EER algorithm in two and three dimensions. 
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d 










' mt,t 




Zint 


2.227 ±0.002 


2.182 ±0.002 


2.198±0.002 


1.673 ±0.001 


2 


B 


0.1100 ±0.0010 


0.0956 ±0.0008 


0.0856 ±0.0008 


0.392 ±0.002 






0.660 


1.04 


0.054 


1470 




Zint 


2.088 ±0.002 


2.052 ±0.002 


2.062 ±0.002 


1.681 ±0.001 


3 


B 


0.1820 ±0.0020 


0.1386±0.0012 


0.1338 ±0.0010 


0.392 ±0.002 






0.804 


0.625 


0.191 


851 



Table 5: Dynamic exponent Zint for the EER algorithm in two and three dimensions, obtained 
by fitting rint = BN^'^'K The number of degrees of freedom of the fit is 1. 



d 






Rl 




£ 


2 


b 


2.21 ±0.02 


2.20 ±0.02 


2.22 ±0.02 


2.14±0.02 




a 


0.0010±0.0010 


0.0025 ±0.0020 


0.002 ±0.002 


0.285 ±0.025 


3 


b 


2.085 ±0.015 


2.068±0.015 


2.085 ±0.015 


2.030 ±0.025 




a 


0.0012 ±0.0010 


0.0025 ±0.0025 


0.0015 ±0.0015 


0.210±0.020 



Table 6: Dynamic exponents a and b for the EER algorithm in two and three dimensions. 



The dynamic critical exponents can also be determined by analyzing the scaling behavior 
of the autocorrelation function, see Eq. (p2D, i.e. by determining a and b so that N'^'^pAAit) 
is a universal function of tN~^, independent of A^. In Table ^ we report the values of a and 
b for which a collapse of the autocorrelation functions is observed. In Fig. ^ we show the 
corresponding plots for i?^ and £ in two dimensions. 

For the energy we observe a very good collapse, while for the radii the scaling behavior 
deteriorates as tN~^ increases. Since the vertical scale changes by several orders of mag- 
nitude, the deviations are somewhat difficult to observe and this makes difficult to set the 
errors on a and b. For this purpose, we found more useful to consider, instead of N"-^ pAAi^)-, 
the quantity 

that is also a universal function of tN~^ in the scaling limit and that scales as t^^^^aN^^ for 
tN-^ oo. 

In Fig. 1^ we report the quantity rscai,yi(^; N) for i?^ and £. The energy shows a very good 
scaling behavior while the scaling is quite poor for i?^, although it improves as increases: 
the data for = 100 and = 300 overlap up to tN'^ ^ 0.05, while the data for A^ = 300 
and A^ = 1000 overlap up to tA^"^ ^ 0.09. 

In three dimensions, all observables show a very good scaling behavior, as it can be seen 
from Fig. ^: In all cases the results for the three values of A^ fall onto a single curve. 

From the results of Table || we can compute the exponents Zexp and Zint, cf. Eqs. (^) 
and (0), and compare them with the previous results. In three dimensions, the estimates 
of Zexp obtained from the scaling analysis are in perfect agreement with those of Table 
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□ N = 100 
A N = 300 
o N = 1000 



0.00 0.10 0.20 0.30 0.40 0.50 
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Figure 6: Dynamic scaling analysis for the EER algorithm in two dimensions: plots of 
InlN'^^PAAit)] vs. tN-K Left frame: A ^ R^, a ^ 0.001, b = 2.21; right frame: A = S, 
a = 0.285, 6 = 2.14. 




Figure 7: Dynamic scaling analysis for the EER algorithm in two dimensions: plots of 
Tscai,A(^; N) vs. ^iV-^ Left frame: A ^ R^, a ^ 0.001, b = 2.21; right frame: A ^ £, 
a = 0.285, 6 = 2.14. 
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Figure 8: Dynamic scaling analysis for the EER algorithm in three dimensions: plots of 
rscai,A(t;iV) vs. tN-\ Left frame: A = R^, a = 0.0012, h = 2.085; right frame: A = £, 
a = 0.210, 6 = 2.03. 



^ In two dimensions instead, only Zexp,£ is compatible with the results of Table ^ The 
estimates of z^xp for the radii obtained from the scaling analysis are significantly lower than 
those obtained from fitting the autocorrelation times. The origin of the discrepancy can 
be understood from Fig. |^. The exponential time is determined by the large-t behavior of 
PAA^t) and in practice by the behavior in the region in which t ^ 3rexp,_R2, see Fig. ||, which 
corresponds approximately to tN^^ ^ 0.44. But in this region there is no scaling, and thus 
the corresponding T^xp.R^ do not scale as N'^^^p-R^ . For this reason, we believe the estimates of 
Table ^ to be grossly in error. Note also that the estimate we obtain in the scaling analysis, 
Zexp,R2 ~ 2.21, is compatible with Zexp,R2 ~ Zint,R2, a relation that is expected to be true, 
since the radii are strongly coupled to the slowest modes of the dynamics. 

Then, we can compute ^mt.A- For the radii, we always have a 0, so that -Zint.^a ~ -2exp,i?2 
as expected. This confirms the results of Table ^ For the energy, we have instead a ^ 0. 
Using Eq. (^4|), we obtain 

;Zint,f = 1.53 ±0.06 d = 2, (29) 

2int,f = 1.60 ±0.05 d = 3. (30) 

These results are in reasonable agreement with those reported in Table ^ if one takes into 
account that, as we discussed, the error on those results is of order 0.04. 

In conclusion, putting the results of the different analyses together, we obtain in two 
dimensions: 

^exp,R2 = 2:int,_R2 = 2.2 ± 0.1, 

2exp,£ = 2.15 ±0.05, 
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-int,£ 



1.60 ±0.10. (31) 



In three dimensions we have: 



^exp,i?2 = ^mt,R2 = 2.07 ± 0.05, 

Ze.p,e = 2.05 ± 0.05, 

Zint,s = 1-65 ± 0.05. (32) 

The errors are such to include the results of all analyses. 

Note that Zi„t,e < ^exp.f for the EER dynamics. This may be understood as follows. The 
energy fluctuations are essentially due to two causes. First, there are fluctuations due to local 
changes of the walk. These fluctuations are fast since are due to local and bilocal moves. 
Then, there are fluctuations due to changes of the global structure of the walk. Indeed, 
there are contributions to the energy that are due to groups of monomers that are far apart 
along the walk but that are near in position space. Such contributions vary slowly, typically 
as rexp,_R2, and are the origin of the fact that Zexp,£ ~ ^exp,/?^. However, these contributions 
are very small, and thus give rise to tiny fluctuations that are negligible when considering 
integrated quantities. Therefore, Zi^t^s < Zcxp,£- 

All results we have discussed up to now refer to simulations with p = 0.5 and with 
the first version of the reptation move. These choices are of no relevance for the critical 
exponents, but they have of course a strong influence on the amplitudes. We have thus tried 
to see the changes in the dynamics due to a variation of p and to a change of the reptation 
move. To check the role of p, we have performed simulations with p = 0.1 and p = 0.9 with 
walks of length N = 100 in two dimensions. We find: for p = 0.9, rint,ij2 = 2184 ± 32 and 
Ti^t,£ = 1200 ± 13; for p = 0.1, T^^t^R2 = 8740 ± 260 and Tint,e = 1181 ± 13. This should be 
compared with the results of Table ^, Tint,/?! ~ 3110 and Tiat,e ~ 862. Thus, by increasing p, 
there is a significant speed up of the dynamics of the radii — this should be expected, since 
reptation moves are essentially the only ones that change the position of the endpoint and 
are thus those that control the slowest modes of the dynamics. Thus, for noninteracting 
SAWs, for which the energy is not an interesting observable, p close to one — but not p = 1, 
otherwise ergodicity is lost — is a good choice for a fast dynamics. On the other hand, the 
dynamics of the energy becomes slower both for p = 0.1 and for p = 0.9. The fact that 
Tint,£ is larger for p = 0.9 is easy to understand. Indeed, by increasing p we decrease the 
probability of performing LO and B22 moves that should be the most important ones for the 
energy. However, it is clear that also reptation moves are relevant for the energy, since for 
p = 0.1 T[at,£ is also larger. Apparently, for small p the relevant quantity for the dynamics 
of the energy is the number of successful moves, irrespective of the type. Indeed, using the 
results of App. |A.lj and |A.2| we find Ti^t,£ ~ 600 successful iterations both for p = 0.5 and 



p = 0.1. 

We have also tested the second version of the reptation algorithm, see App. |A.2| . For 



random walks, this implementation gives r ~ compared to r ~ N"^ of the first version. For 
the SAW, the two versions are expected to have the same critical exponents, but the second 
one should be faster. We have performed a simulation with p = 0.5 in two dimensions, finding 
Tint,_R2 = 510 ± 4 and rint,£- = 299.8 ± 1.6. These estimates are sensibly smaller than those 
reported in Table H. For the dynamics is faster by a factor of 6, and for £^ by a factor of 3. 
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N 


iter 




'^exp,i?2 


'^exp,i?^ 


100 


2.4- 10" 


37600 ± 160 


38000 ±160 


38620 ±100 


300 


8.8-10" 


918000 ±8000 


958000 ±8000 


972000 ±8000 


700 


1.9- 10^2 


11475000 ±315000 


11610000 ±315000 


11250000 ±225000 



Table 7: Exponential autocorrelation times for the KER algorithm in two and three dimen- 
sions, "iter" is the number of iterations. 





■7"cxp,i?2 


■^cxpjiJ^ 


■7"cxp,i?2„ 




2.92 ±0.01 


2.94 ±0.01 


2.93 ±0.01 


B 


0.054 ±0.002 


0.050 ±0.002 


0.054 ±0.002 




3.77 


0.034 


2.46 



Table 8: Dynamic exponent z^^p for the KER algorithm in two and three dimensions, ob- 
tained by fitting Tgxp = BN^"^'^ . The number of degrees of freedom of the fit is 1. 



Clearly, the second implementation is the most efficient one and all our simulations should 
have used it. Unfortunately, we thought of this second version only when all simulations 
were completed. 

6 The KER dynamics in two dimensions 

The second dynamics we have analyzed is the KER algorithm in which end-end reptation 
moves BEE are replaced by kink-end reptation moves BKE. It turns out that this dynamics 
is much slower than the EER one, and thus we have limited our analysis to shorter walks, 
= 100, 300, 700, to two dimensions, and to noninteracting SAWs, i.e. /5 = 0. We set 
p = 0.5. Again, we measured three radii and the energy. The static results agree with those 
obtained by using the EER dynamics and discussed in the preceding section and with the 
results of Ref. jSTj. 

As in the preceding section, we have first determined the exponential autocorrelation 
times by studying the large-time behavior of the effective exponents fexp,A(t', s), see Sec. |[ 
The results for the radii are reported in Table |^. We have not been able to determine the 
exponential autocorrelation time for the energy S. Indeed, for the values of t for which pes{t) 
is not zero within statistical errors, pss{i) has a power-law behavior, i.e. pse{i) ~ ^ with 
a ~ 1-1.3. This can clearly be seen from Fig. |^ where we report pes{i) for N = 300. 

The results for the exponential autocorrelation times are fitted with the Ansatz (P^D, 
obtaining the results reported in Table ^ As before, we cannot perform a systematic analysis 
of the scaling corrections. However, it is important to notice that the estimates for the three 
radii agree within error bars. This confirms the correctness of the quoted error bars and 
gives the final estimate 

-2exp,i?2 = 2.93 ± 0.02. (33) 
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This estimate is significantly higher than that for the EER dynamics. The origin of such a 
large difference is unclear, since it is difficult to see any difference between BEE and BKE 
moves. Naively, one would have expected a BKE move to be equivalent to two BEE moves 
together with a B22 move. Since all moves have a finite probability of success as N oo, 
see Appendix, one would have expected a difference by a constant factor, and thus the same 
critical exponents. Such a naive expectation is not true, since the exponents are clearly 
different. 

We have then determined the integrated autocorrelation times. For the radii we have 
used the self-consistent windowing method of Sec. using c = 15. The results are reported 
in Table ^. Since roxp,R2 ~ Srint,/?^, the choice c = 15 should be enough to avoid systematic 
errors. Instead, the autocorrelation function of the energy decreases rapidly and it has a very 
long tail. In this case, we have chosen a much larger valus of c, c = 200, obtaining the results 
that are reported in Table ^ as rint,£-(n.t.). In spite of c being such a large number, the cutoff 
value M, cf. Eq. (JT^), is still well within the region in which the function decays as a power 
law. For instance, for = 300, M = 670000 (InM ^ 13.4) and p{M) ^ 3.4 ■ 10"^, see Fig. 
^ This should be expected since the cutoff M satisfies M < rexp,_R2 ~ 10^. More precisely, 
M ^ 3.3rexp,ij2, M ^ 0.7rexp,R2, and M ^ 0.14rexp,R2 approximately for N = 100,300,700. 
Therefore, we expect a sizable contribution from the tail of the autocorrelation function, at 
least for = 300, 700. To take it into account, we use Eq. (pOD, where for rexp,£: we take the 
average rexp,_R2 of the radii. For A^ = 300 the correction is of 2.5% (2.4% if we approximate 
Peeif) ~ Bt~^-^ as obtained from the fit) and for A^ = 700 of 20.3% (18.1% if we approximate 
Pseif) ~ Bt~^-^ as obtained from the fit). The results are reported in Table |] as rint,£-(w.t.). 




168000 336000 504000 672000 ^ ,, , 

t Log(t) 



Figure 9: Autocorrelation function pssit) for the KER algorithm in two dimensions. Here 
A^ = 300. The straight line in the right frame corresponds to Bt~°', B = 5100, a = 1.29. 
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N 












100 


15632 ±30 


20974 ±48 


27528 ±72 


847.5 ±0.4 


847.5 ±0.4 


300 


341030 ±1640 


462700 ±2600 


648500 ±4300 


3402.6 ±1.6 


3484 ±4 


700 


3906500 ±43800 


5325000 ±70000 


7942000 ±127000 


9269.0±5.1 


10948 ±90 


Zint 


2.820 ±0.005 


2.830 ±0.004 


2.890 ±0.006 


1.2330 ±0.0004 


1.289 ±0.0011 


B 


0.0360 ±0.0007 


0.0464±0.0010 


0.0458 ±0.0012 


2.920 ±0.006 


0.808 ±0.005 




18.5 


11.8 


11.8 


4197 


41.7 



Table 9: KER algorithm in two dimensions: integrated autocorrelation times and dynamic 
critical exponents z^-at-, obtained by fitting Tint = BN^^'^^. The number of degrees of freedom 
of the fit is 1. For the energy £, we report two estimates: (n.t.) is obtained by using the 
self-consistent windowing method with c = 200, while (w.t.) is the result obtained including 
the tail contribution a la Li et al. (Ref. ||57|). 









£ 


b 2.89±0.03 


2.92±0.02 


2.92±0.02 


[2.93 ±0.04] 


a 0.015±0.015 


0.020±0.015 


0.015±0.010 


0.670±0.015 



Table 10: Dynamic exponents a and b for the KER algorithm in two dimensions. For the 
energy, we have fixed b = 2.93 ± 0.04 and determined the corresponding a. 



The reported error is obtained by summing 20% of the contribution of the tail to the original 
error. This is an "ad hoc" prescription which can be shown to work reasonably in the exactly 
soluble case of the pivot algorithm for the random walk. 

The results for the autocorrelation times have been fitted with the Ansatz (^) in order to 
obtain Zjnt- For the radii, the quality of the fits is quite poor, with a of approximately 10- 
20. Moreover, the estimates do not agree within error bars. There are therefore corrections 
to scaling larger than the very tiny statistical errors. By requiring Ziat,R^ to coincide for these 
three quantities, we obtain finally 

2int,ij2 = 2.85 ± 0.06, (34) 

that includes all estimates and is compatible with the expectation 2;int,R2 = Zexp,R2 ■ 

Fits of Ti^t,£ are characterized by a very large ^ind give 2int,£: ~ 1.2-1.3, much smaller 
than Z[nt,R^- As in the EER algorithm, the dynamics of the energy is much faster than that 
of the radii. Note that -Zint.f is also significantly lower that the corresponding exponent for 
the EER algorithm. Again, it is quite difficult to understand intuitively why this happens. 

The results reported above are confirmed by a scaling analysis using Eq. (0). The 
results for the radii are reported in Table |l^. The exponent a is compatible with zero and 
^exp,_R2 = b = 2.91 ± 0.03, in agreement with the estimate (p^). 

Since p£e{t) ~ Bt'" for the values of t we can investigate, we cannot determine a and 
b independently. Thus, we have assumed Zexp,e = Zexp,R2 and then used b = 2.93 ± 0.04, 
cf. Eq. ( p3D — to be conservative, we have doubled the error. Correspondingly, we obtain 
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TV 


iter 






'^exp,i?^ 




100 


5.84- IQi" 


5120 ±70 


5220 ±60 


5170 ±40 


4856 ±44 


800 


6.28- 10" 


609600 ± 15000 


606000 ±10000 


620000 ±8000 


594200 ±9000 


1600 


1.81 • 10^2 


3165000 ±22000 


2840000 ±20000 


3020000 ±20000 


2936000 ±24000 


3200 


1.12 • 1013 


15838000 ±76000 


13880000 ±50000 


14990000 ±50000 


14656000 ±50000 



Table 11: Exponential autocorrelation times for the EER algorithm in two dimensions at 
the 6 point, "iter" is the number of iterations. 



a = 0.670 ± 0.015 and Zint,£ = 0.97 ± 0.05, which is somewhat lower than the estimates of 
Table |. One may think that this is due to our assumption Zexp,s = ^exp,_R2, while there is 
some evidence from the analysis of the EER dynamics that Zexp,£ < ^exp.R^- However, this 
does not explain the difference, since if Zexp,£ decreases, also Zint,^ decreases. In order to 
obtain Zi^t,£ = 1-3, one should take Zexp,£ = b = 3.35, which is much too large. Therefore, 
the difference should be taken as an indication of the scaling corrections. 

In conclusion, the KER dynamics has a different critical behavior with respect to the 
EER dynamics. For the critical exponents we have 

^exp,ij2 = Z;^t,R^ = 2.90 ± 0.05, (35) 
Zint,^ = 1.0 ± 0.3. (36) 



7 The EER dynamics at the point in two dimensions 

Bilocal algorithms are of interest for applications in constrained geometries and in the pres- 
ence of strong interactions where nonlocal algorithms are inefficient. 

In this section we study the dynamic behavior of the EER algorithm at the 6 point in 



two dimensions, by setting (3 = /Sg = 0.665 — we have used the estimate of Ref. |^, see 



Table |I|. Here, we have studied longer walks than before, N = 100, 800, 1600, 3200, with 



large statistics, see Table ITT 



We have performed the same analyses we have presented in the preceding sections. First, 
we have determined the exponential autocorrelation times. For all observables, the effective 
exponent fexp,yi(^; s) becomes independent of tfoit^ 2-2.5 r^^p^A, allowing a reliable estimate 



of the exponential autocorrelation times. The results are reported in Table |TT. 

Then, we determined the exponent Zexp^A by fitting the exponential autocorrelation times 
to the Ansatz (^Sf) . The results are reported in Table |12|. Clearly, the statistical errors are 
too small. Indeed, we expect Zexp,A to be the same for all observables — including the energy, 
that should be strongly coupled to the slowest modes at the 9 point — and this does not 
happen with the quoted error bars. By direct comparison of all estimates, we obtain the 
more conservative result 

Zexp,£ = Zexp,R2 = 2.30 ± 0.03, (37) 

where the error is such to include all estimates. 
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exp,i?^ "^exp^Rl '^exp,fl^ "^exp.g 



2cxp 


2.318±0.004 


2.275 ±0.003 


2.300 ±0.004 


sr 1 

2.313 ±0.004 


B 


0.126±0.006 


0.148 ±0.004 


0.128 ±0.005 


0.116±0.006 




2.01 


3.64 


1.24 


0.78 



Table 12: Dynamic exponent z^^^ for the EER algorithm in two dimensions at the 6 point, 
obtained by fitting Tcxp = -BA^^^^p. The number of degrees of freedom of the fit is 2. 



N 


^int,i?2 


''"int,i?2 






100 


4520 ± 10 


3166.2 ±5.8 


3082.8 ±5.4 


2504.4 ±4.0 


800 


559500 ±4100 


362300 ±2100 


357900 ±2100 


222800 ±1000 


1600 


2890400 ±28300 


1779500 ±13700 


1778000 ±13600 


1029500 ±6000 


3200 


13953000 ±120000 


8630000 ±59000 


8698000 ±59000 


4185000 ±20000 




2.321 ±0.002 


2.282 ±0.002 


2.291 ±0.002 


2.1510±0.0011 


B 


0.1020±0.0013 


0.0864 ±0.0008 


0.0800 ±0.0010 


0.1260±0.0008 




9.0 


1.22 


3.73 


149 



Table 13: EER algorithm in two dimensions at the 6 point: integrated autocorrelation times 
and dynamic critical exponents ^int, obtained by fitting Tint = BN^'^^'K The number of degrees 
of freedom of the fit is 2. 



We have analogously determined the integrated autocorrelation times using the self- 
consistent windowing method with c = 15. The results for rint,A are reported in Table |13 



Notice that in this case the integrated autocorrelation times for the energy are close to those 
of the radii, as it should be expected, since at the 9 point also the energy is a "slow" variable. 
In all cases, t^^^ a ~ 1-3 Tint,^ and thus the choice c = 15 should give a small systematic error 
due to the truncation of the autocorrelation functions. 

The integrated autocorrelation times have been fitted with the Ansatz i?A^^'"% in order 
to compute Zint,^- The results are reported in Table |I3|. In all cases, the purely statistical 
errors we have reported are too small. For the radii, the exponent z^^i should be the same. 



and thus the error is at least a factor of ten larger. Comparing the estimates of Table |13 
we arrive at the final result 

2.30 ±0.03, (38) 



int,ii;2 



that, by comparing with Eq. 



gives 2;int,ij2 = 2exp,i?2, as expected. 



The result for E is somewhat lower, but the very large indicates that corrections to 
scaling are significant. In order to see if there is a systematic trend we have determined an 
effective ^^mt/, by computing 



In 



rint,g(A^i 



N2 



(39) 



We obtain 



(100,800) 



2.158 ±0.003, 



(40) 
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Zint,£: (800, 1600) 
Zint,£(1600,3200) 



2.208 ±0.015 
2.023 ±0.015 



(41) 
(42) 



It is difficult to observe a systematic trend, but in any case a systematic increase towards 
2.30 seems to be excluded. On the contrary, the data seem to indicate that Zij^t,£ decreases 
below the value of Table |TB|. Thus, also at the 9 point we have Zmt.s < Ze^p^e, although the 
difference is much smaller than in the case (3 = 0. 



These results are con&med by a scaling analysis. In Fig. ITO we report the scaling variable 



Tscai^Ri^it; N). We observe a very good scaling behavior and correspondingly we are able to 
obtain quite reliable estimates of the critical exponents a and b. The same good behavior 
is observed for all observables. The estimates of a and b are reported in Table 11^. For all 



observables, b is compatible with the estimates of Table [T^, conffiming the estimate ([371) . 
For the radii, a = 0, in agreement with Table [T^ For the energy, a is clearly nonvanishing, 
conffiming that -Zmt.f < -Zexp,£:- Using Eq. (p^), we have ^int,£- = 2.15 ± 0.03, which is in 
agreement with the previous results. 




Figure 10: Dynamic scaling analysis for the EER algorithm in two dimensions at the 9 point: 
plots of rscai,i?2^(t; N) vs. tN~K Here a = 0.0015, b = 2.31. 



22 







ry-2 


£ 


b 2.315±0.015 


2.30±0.02 


2.31 ±0.02 


2.31 ±0.02 


a 0.0015 ±0.0015 


0.002 ±0.002 


0.0015 ±0.0015 


0.07 ±0.01 



Table 14: Dynamic exponents a and b for the EER algorithm in two dimensions at the 6 
point. 



8 Conclusions 

The simulations we have presented show that the reptation dynamics is quite successful, even 
at the 6 point. The values of z we have found are only marginally higher than 2, which is the 
best possible behavior for a dynamics that involves local and bilocal moves. Of course, for 
a practical implementation one may want to explore several variants that, although do not 
change the critical behavior, may still speed up the dynamics by a constant (large) factor. 
First of all, in practical implementations it is important to use the second version of the 
reptation dynamics (see App. |A.2|) . Second, the B22 moves are quite rarely performed and 
in any case much less than the kink-end/end- kink moves. For instance, in two dimensions at 
/3 = 0, B22 (resp. BKE) moves are performed with probability 0.08 (resp. 0.13). Moreover, 
BKE moves appear to be quite successful in speeding up the dynamics of the energy, that is 
one of the slow variables in the presence of interactions. Therefore, in the compact regime an 
efficient dynamics can be obtained by mixing together: (i) the reptation move; (ii) the BKE 
move; (iii) purely local moves LO and LI. A purely local algorithm that leaves the correct 
measure invariant can be obtained from that described in App. |A.1| by setting p{0) = 1/2 
in all dimensions and p(22) = 0. In this case, it is convenient to include also crankshaft 
moves, see Sec. 4.1 of Ref. |T^. Such an implementation of the EER algorithm should be 
the method of choice for fixed simulations in the compact regime. 



A The basic moves 

In this appendix we introduce the basic moves that we use in our simulation: 

(i) The kink- kink local/bilocal move; 

(ii) The reptation move; 

(iii) The kink-end/end- kink reptation move. 



In Ref. [16] it was shown that moves (iii) are enough to obtain an ergodic algorithm. In two 
dimensions one can limit oneself to consider only moves (i), but this algorithm is inefficient 
because of the slow motion of the endpoint. Reptation moves are never ergodic because of 
the possibihty that the endpoints get trapped. 
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(a) 



(b) 



(c) 



(d) 



Figure 11: Configurations of tliree consecutive links: (a) configuration of type I; (b) confi^ 
uration of type L; (c) configuration of type S; (d) configuration of type U. 



A.l Kink-kink local/bilocal move 



In tliis section we define tlie kink-kink local/bilocal move |jT6[. In order to describe the 
algorithm it is important to classify the possible configurations of three successive links (see 



Fig. |n 



1. the bonds have the same direction (I configuration); 

2. two consecutive bonds have the same direction, while the third one is perpendicular to 
them (L configuration); 

3. the first and the third bond are perpendicular to the second one, and they are either 
parallel or perpendicular to each other (S configuration); 

4. the first and the third bond are perpendicular to the second one, and they are antipar- 
allel to each other (U configuration). 

An iteration works as follows: 

• Step 1. Choose a random site i of the current walk a;, < i < iV. li i = N, propose 
an LI move and go to step 5. 

• Step 2. Determine the configuration of the subwalk uj[i — l,i + 2]. If i = N — 1, we 
imagine adding a link Alj[N) parallel to Auj{N — 1), so that the possible configurations 
are of type L and I. Analogously, if i = 0, we imagine adding a link Aa;(— 1) parallel 
to Acj(O). 

• Step 3. Draw a random number r, uniformly distributed in [0,1]. Depending on the 
configuration of uj[i — 1, i + 2], do the following: 
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1. I: If r > (2(i — 2)p(22), perform a null transition and the iteration ends. Otherwise, 
go the next step. 

2. L: If r > {2d — 3)p(22) +p(0), perform a null transition and the iteration ends. If 
{2d — 3)p(22) < r < {2d — 3)p(22) +p(0), propose an LO move and go to step 5. 
Otherwise, go to the next step. 

3. S: If r > {2d — A)p{22) + 2p{0) perform a null transition and the iteration ends. 
If {2d — 4)j9(22) < r < {2d — 4)p(22) + 2p(0) propose an LO move: there are two 
possibilities which are chosen amongst randomly; then go to step 5. Otherwise, 
go to the next step. 

4. U: Go to the next step. 

• Step 4. Choose a second integer j uniformly in the disjoint intervals, — 1 < j < A^, 
j ^ i — + If j = —1,N make a null transition and the iteration ends. Otherwise, 
depending on the configuration of uj[i — l,i + 2], do the following: 

— u[i — l,i + 2] is of type I, S, L: if j = or j = iV — 1, or if u![j — 1, j + 2] is not 
of type U perform a null transition and the iteration ends. Otherwise, propose a 
B22 move, cutting the kink uj[j — 1, j + 2] and adding it to uj[i,i + 1] in one of the 
possible directions ||6^. Then, go to the next step. 

— u[i — l,i + 2] is of type U: according to the configuration of uj[j — l,j + 2] (if 
j = 0, — 1 imagine adding links as before) do the following: 

1. uj[j — 1, j + 2] is of type I: If r < {2d — 2)p{22) (note that the random number r 
appearing here is the same used in Step 3.), propose a B22 move: cut the kink 
uj[i — l,i + 2] and add it on top of uj[j,j + 1] in a possible random direction, 
and then go to step 5. Otherwise, perform a null transition and the iteration 
ends. 

2. a;[j — 1, j + 2] is of type L: If r < {2d — 3)p(22), propose a B22 move: cut 
the kink u![i — l,i + 2] and add it on top of uj[j,j + 1] in a possible random 
direction, and then go to step 5. Otherwise, perform a null transition and 
the iteration ends. 

3. uj[j — 1, j + 2] is of type S: If r < {2d — 4)p(22) propose a B22 move: cut 
the kink u}[i — l,i + 2] and add it on top of uj[j,j + 1] in a possible random 
direction, and then go to step 5. Otherwise, perform a null transition and 
the iteration ends. 

4. uj[j — 1, j + 2] is of type U: If r < {2d — 3)p(22), propose a B22 move: cut 
the kink uj[i — l,i + 2] and add it on top of uj[j,j + 1] in a possible random 
direction, and then go to step 5; if {2d — 3)p{22) < r < 2{2d — 3)p{22), propose 
a B22 move: cut the kink uj[j — 1, j + 2] and add it on top of uj[i, i + 1] in a 
possible random direction, and then go to step 5. Otherwise, perform a null 
transition and the iteration ends. 

• Step 5. Check for self-avoidance. If the proposed new walk is self-avoiding keep it, 
otherwise perform a null transition. 
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d 


N 




LO, LI 


B22 




100 


Ppr 


0.496 


0.0952 






PSA 


0.751 


0.852 




300 


Ppr 


0.499 


0.0963 


2 




PSA 


0.751 


0.851 




700 


Ppr 


0.500 


0.0966 






PSA 


0.751 


0.851 




1000 


Ppr 


0.500 


0.0967 






PSA 


0.751 


0.851 




100 


Ppr 


0.443 


0.0850 






PSA 


0.800 


0.882 


3 


300 


Ppr 


0.446 


0.0864 






PSA 


0.798 


0.878 




1000 


Ppr 


0.447 


0.0869 






PSA 


0.798 


0.877 



Table 15: Proposal probability Ppj. and probability psA that the proposed walk is self-avoiding 
for local (LO and LI) and bilocal (B22) moves. We consider noninteracting SAWs in two 
and three dimensions. 



• Step 6. Compute the difference in energy between the old and the new walk and 
perform a Metropolis test. 

The algorithm we have presented depends on the probabilities p{0) and J9(22), that are 
the probabilities of LO and B22 moves respectively. As discussed in Ref. |TB|, the fastest 
dynamics is obtained by setting 

In two dimensions p(0) = p{22) = 1/2, while in three dimensions p(0) = l/3andp(22) = 1/6. 

It is interesting to compute the probability of a successful move. For noninteracting 
SAWs such a probability is the product of two terms: the probability that a given move 
is proposed and the probability psA that the proposed walk is self-avoiding. At the 6'-point 
one must additionally multiply by the probability that the Metropolis test is successful. 
Numerical estimates of these probabilities are reported in Tables ^ and |16[ Note that they 
have a very weak N dependence and clearly approach a constant value as goes to infinity. 

The probability can be easily computed by using the probability of occurrence of the 
four configurations I, L, U, and S defined in Fig. 0. Indeed, 

Ppr,B22 = 2p(22MU)[(2rf-2MI) + (2rf-3)(j9(L)+p(U)) + (2rf-4)p(S)], (45) 
Ppr,L = p(0)[p(L) + 2p(S)]. (46) 



Using the numerical results of Table |T7|, we obtain pp^ = 0.097, 0.087, 0.164 for B22 moves 
and {d = 2, P = 0), {d = 3, P = 0), and {d = 3, P = Pe) respectively. For local moves we 
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N 




LO, LI 


B22 




Ppr 


0.463 


0.158 


100 


PSA 


0.429 


0.525 




PMct 


0.701 


0.610 




Ppr 


0.462 


0.162 


800 


PSA 


0.394 


0.416 




PMct 


0.711 


0.561 




Ppr 


0.462 


0.163 


1600 


PSA 


0.389 


0.399 




PMct 


0.714 


0.556 




Ppr 


0.462 


0.163 


3200 


PSA 


0.385 


0.388 




PMct 


0.716 


0.553 



Table 16: Proposal probability ppr, probability psA that the proposed walk is self-avoiding, 
and probability pMet that the proposed self- avoiding walk is accepted in the Metropolis test. 
For local (LO and LI) and bilocal (B22) moves. We consider SAWs in two dimensions at the 
9 point. 





13 = 0, d = 2 




(3 = 0, d = 3 


1 


0.151 


0.128 


0.051 


L 


0.480 


0.455 


0.354 


U 


0.109 


0.183 


0.102 


S 


0.260 


0.234 


0.493 



Table 17: Probabilities for ^ oo of the occurence of the four configurations I, L, U, and 
S defined in Fig. |lT[ Results for noninteracting SAWs in two and three dimensions and for 
SAWs at the 6 point in two dimensions. 



obtain correspondingly ppj- = 0.500, 0.462, 0.447. These results are in good agreement with 



the numerical ones of Tables O and 16 



Using the above presented results, we can compute the probability of a successful move. 
They are reported in Table |18[ Note that at the 6 point the probability of a null transition 
is quite large and in particular B22 moves are quite rarely performed. 



A. 2 Reptation move 

There are two different implementation of the reptation (or slithering-snake) move. The first 
one, which satisfies detailed balance, works as follows [Version 1]: 

• Step 1. With probability 1/2 delete uj[N — 1, N] and add a new link at the beginning 
of the walk; otherwise, delete ^^[0, 1] and add a new link at the end of the walk. 

• Step 2. Check if the new walk is self- avoiding. If it is keep it, otherwise perform a null 
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d 


/? 


L0,L1 


B22 


null 


2 





0.38 


0.08 


0.54 


3 





0.36 


0.08 


0.56 


2 


Pe 


0.13 


0.03 


0.84 



Table 18: Probability of the different moves for different (3 and d. 



d 


N = 100 


N ^ 300 


N = 1000 


2 


0.882 


0.880 


0.880 


3 


0.938 


0.937 


0.937 



Table 19: Probability psA thai the proposed walk is self-avoiding for the reptation move in 
two and three dimensions. We consider noninteracting SAWs. 



transition. 

• Step 3. Compute the difference in energy between the old and the new walk and 
perform a Metropolis test. 

A second version uses an additional flag which specifies which of uj{0) and uj{N) is the 
"active" endpoint. It works as follows [Version 2]: 

• Step 1. Delete one bond at the "active" endpoint and append a new one at the opposite 
end of the walk. 

• Step 2. If the new walk is self-avoiding keep it, otherwise stay with the old walk, and 
change the flag, switching the active endpoint. 

• Step 3. Compute the difference in energy between the old and the new walk and 
perform a Metropolis test. 

This algorithm |^ does not satisfy detailed balance, but it satisfies the stationarity condition 
generating the correct probability distribution. 





N = 100 


N ^ 800 


N = 1600 


N = 3200 


PSA 


0.643 


0.566 


0.551 


0.540 


PMct 


0.697 


0.663 


0.658 


0.654 



Table 20: Probability psA that the proposed walk is self- avoiding and probability pMet that 
the proposed self-avoiding walk is accepted in the Metropolis test. For reptation moves. We 
consider SAWs in two dimensions at the 9 point. 
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It is interesting to compute the probability of success of a reptation move. In the absence 
of interactions it is simply given by the probability that the proposed walk is self-avoiding. 
Such a probability is reported in Table |19[ The reptation move is quite successful, being 
accepted with high probability in both two and three dimensions. 

At the 9 point, we must also consider the probability that the proposed walk passes the 
Metropolis test. Numerical results are reported in Table |20|. Since the walk is more compact, 
PsA is lower than in the noninteracting case, although still quite large. Multiplying the two 
probabilities we see that the reptation move is accepted in 35% of the cases. Note that this 
probability is larger than the probability of a local or bilocal B22 move, see Table 



A. 3 Kink- end /end- kink move 

The kink-end/end-kink move uses BKE moves (see Fig. |^). It consists of the following steps: 

• Step 1. Choose a random site i of the current walk with < z < — 2. 

• Step 2. Propose an end-kink move with probability {2d — 2)p or a kink-end move with 
probability {2d — the first case delete the last two bonds of the walk and insert 
a kink on the bond Auj{i) in one of the {2d — 2) possible orientations. In the second 
case, if i 7^ and uj[i — l,i + 2] is a kink, remove it and attach two bonds at the end 
of the walk in one of the {2d — 1)^ possible ways. Otherwise, perform a null transition 
and the iteration ends. 

• Step 3. Check if the proposed walk is self-avoiding. If it is keep it, otherwise make a 
null transition. 

• Step 4. Compute the difference in energy between the old and the new walk and 
perform a Metropolis test. 

The constant p is given by [|l^ 

^ " {2d - ly + {2d - 2)' ^^^^ 

We obtain p =l/llin(i = 2, andp = l/29in(i = 3. A slightly more efficient implementation 
is discussed in Ref. |]16| . 

We have computed numerically the probability that a kink-end or an end-kink move is 
accepted. We find 0.140, 0.138, 0.138 for A^ = 100,300,700 respectively. The probability of 
a null transition is therefore quite large, much larger than for a kink-kink bilocal move. Note 
however, that a kink-end/end-kink move is performed more often than a B22 move and thus 
this type of moves should be slightly more efficient in updating the part of the walk that is 
far from the endpoints than B22 moves. 
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